CC para resíduos βAR(1)

cache_dados <- function(chave, funcao_geradora) {
  cache <- file.exists(chave)
  
  if (!cache) {
    resultado <- funcao_geradora()
    saveRDS(resultado, chave)
  } else {
    resultado <- readRDS(chave)
  }
  
  resultado
}

Parâmetros do modelo

nH0 <- 100
nH1 <- 200
phi_parametro <- 0.2

coeficientes <- function(phi) {
  # βARMA: the model from Rocha and Cribari-Neto (2009, 2017)
  #        is obtained by setting coefs$d = 0 and d = FALSE and error.scale = 1 (predictive scale)
  list(alpha = 0, beta = NULL, phi = phi, theta = NULL, d = 0, nu = 20)
}

barma.sim <- function(n, phi, seed, y.start = NULL){
  BARFIMA.sim(
    n = n,
    coefs = coeficientes(phi),
    y.start = y.start,
    error.scale = 1,
    complete = F,
    seed = seed
  )
}

barma.phi_estimado <- function(yt, alpha = 0, nu = 20, phi = 0.1){
  BARFIMA.fit(
    yt = yt,
    start = list(alpha = alpha, nu = nu, phi = phi),
    p = 1,
    d = FALSE,
    error.scale = 1,
    report = F
  )$coefficients["phi"][[1]]
}

barma.residuos <- function(yt, phi_estimado){
  BARFIMA.extract(
    yt = yt,
    coefs = coeficientes(phi_estimado),
    llk = F,
    info = F,
    error.scale = 1
  )$error
}

# Função para forçar a execução de `BARFIMA.extract` até que ela funcione
so_quero_que_funcione <- function(func, debug = T, max = 1000, ...) {
  FINALMENTE <- F
  contador <- 0
  
  while (!FINALMENTE) {
    contador <- contador + 1
    if (debug) print(paste("Tentativa:", contador))
    
    resultado <- tryCatch({
      func(...)
    }, error = function(err) {
      if (contador >= max) {
        stop("Deu ruim, número máximo de tentativas atingido")
      }
      
      if (any(grepl("\\.btsr\\.extract\\(", deparse(err$trace$call)))) {
        # Ignora erro de extração de resíduos
        return(NULL)
      }
      
      stop(err)
    })
    
    if (!is.null(resultado)) {
      FINALMENTE <- T
    }
  }
  
  return(resultado)
}

Carta EWMA

ewma_qcc <- function(amostra_inicial, lambda, amostra_subsequente) {
  ew <- ewma(amostra_inicial, lambda = lambda, nsigmas = 3, plot = F, newdata = amostra_subsequente)
  
  registros <- (nH0 + 1):(nH0 + nH1)
  
  ewma <- as.numeric(ew$y[registros])
  LI <- ew$limits[, 1][registros]
  LS <- ew$limits[, 2][registros]
  fora_de_controle <- ewma < LI | ewma > LS
  total_fora_de_controle <- sum(fora_de_controle, na.rm = TRUE)
  fracao_fora_de_controle <- total_fora_de_controle / length(ewma)
  list(
    ewma = ewma,
    LI = LI,
    LS = LS,
    fora_de_controle = fora_de_controle,
    total_fora_de_controle = total_fora_de_controle,
    fracao_fora_de_controle = fracao_fora_de_controle
  )
}

Simulação de amostras subsequentes

ewma_amostras_subsequentes_gerador <- function() {
  lambda <- 0.2
  
  amostra_controle <- barma.sim(
    n = nH0,
    phi = phi_parametro,
    seed = 1337
  )
  
  ultima_observacao <- amostra_controle[nH0]
  phi_estimado <- barma.phi_estimado(amostra_controle)
  residuos_controle <- barma.residuos(amostra_controle, phi_estimado)
  
  data.frame(
      novo_phi = seq(0, 0.8, by = 0.1)
    ) %>%
    rowwise() %>%
    mutate(
      observacao = list(1:nH1),
      dados = list(barma.sim(
        n = nH1,
        phi = novo_phi,
        seed = 1337,
        y.start = ultima_observacao
      )),
      residuos = list(
        barma.residuos(dados, phi_estimado)
      ),
      qcc = list(ewma_qcc(residuos_controle, lambda, residuos)),
      ewma = list(qcc$ewma),
      LI = list(qcc$LI),
      LS = list(qcc$LS),
      fora_de_controle = list(qcc$fora_de_controle),
      total_fora_de_controle = qcc$total_fora_de_controle,
      fracao_fora_de_controle = qcc$fracao_fora_de_controle
    )
}


ewma_amostras_subsequentes <- so_quero_que_funcione(ewma_amostras_subsequentes_gerador)

## [1] "Tentativa: 1"

ewma_amostras_subsequentes %>%
  group_by(novo_phi) %>%
  summarise(
    mean = mean(fracao_fora_de_controle),
    .groups = "drop"
  )

Cartas

ggplotly(
  ewma_amostras_subsequentes %>%
    unnest(
      cols = c(novo_phi, observacao, ewma, LI, LS)
    ) %>%
    ggplot() +
    geom_line(aes(x = observacao, y = LI, color = "Limite Inferior"), linetype = "dashed") +
    geom_line(aes(x = observacao, y = LS, color = "Limite Superior"), linetype = "dashed") +
    geom_line(aes(x = observacao, y = ewma, color = "EWMA")) +
    geom_point(data = . %>% filter(ewma < LI | ewma > LS),
               aes(x = observacao, y = ewma, color = "Fora de controle"), size = 1.5, shape = 4) +
    labs(x = "Observação", y = "Caracterísitca", color = "Legenda") +
    scale_color_manual(
      values = c(
        "Limite Inferior" = adjustcolor("red", alpha.f = 0.5),
        "Limite Superior" = adjustcolor("red", alpha.f = 0.5),
        "EWMA" = adjustcolor("black", alpha.f = 0.8),
        "Fora de controle" = adjustcolor("red", alpha.f = 0.4)
      )
    ) +
    labs(
      title = "EWMA para resíduos βAR(1), com λ = 0.2",
    ) +
    theme.base +
    theme(
      legend.position = "bottom",
      legend.title = element_blank()
    ) +
    facet_wrap(~novo_phi)
) %>%
  plotly.base

Simulação de Monte Carlo

ewma_monte_carlo_gerador <- function(){
  numero_de_execucoes <- 200
  
  expand.grid(
    k = 1:numero_de_execucoes,
    novo_phi = seq(0.2, 0.6, by = 0.1),
    lambda = seq(0.1, 0.4, by = 0.1)
  ) %>%
  mutate(
    id = row_number()
  ) %>%
  rowwise() %>%
  mutate(
    amostra_controle = list(barma.sim(
      n = nH0,
      phi = phi_parametro,
      seed = id
    )),
    ultima_observacao = amostra_controle[nH0],
    phi_estimado = barma.phi_estimado(amostra_controle),
    residuos_controle = list(
      barma.residuos(amostra_controle, phi_estimado)
    ),
    dados = list(barma.sim(
      n = nH1,
      phi = novo_phi,
      y.start = ultima_observacao,
      seed = id + 1337E4
    )),
    residuos = list(
      barma.residuos(dados, phi_estimado)
    ),
    qcc = list(ewma_qcc(residuos_controle, lambda, residuos)),
    ewma = list(qcc$ewma),
    LI = list(qcc$LI),
    LS = list(qcc$LS),
    fora_de_controle = list(qcc$fora_de_controle),
    total_fora_de_controle = qcc$total_fora_de_controle,
    fracao_fora_de_controle = qcc$fracao_fora_de_controle
  )
}

ewma_monte_carlo <- cache_dados(
  "cache_simulacao_ewma.RData",
  \() so_quero_que_funcione(ewma_monte_carlo_gerador)
)

Resumo

ewma_monte_carlo_resumo <- ewma_monte_carlo %>%
  group_by(lambda, novo_phi) %>%
  summarise(
    mean = mean(fracao_fora_de_controle),
    min = min(fracao_fora_de_controle),
    max = max(fracao_fora_de_controle),
    .groups = "drop"
  )

ewma_monte_carlo_resumo

Gráficos

ggplotly(
  ewma_monte_carlo_resumo %>%
      ggplot(aes(x = novo_phi, y = mean, color = factor(lambda))) +
      geom_line() +
      geom_point() +
      labs(
        x = "Valores de Φ",
        y = "Fração de pontos fora de controle",
        color = "Valores de λ",
        title = "Fração de pontos fora de controle por Φ e λ"
      ) +
      theme.base
) %>%
  plotly.base

Box plot

ggplotly(
  ewma_monte_carlo %>%
    ggplot(aes(x = factor(novo_phi), y = fracao_fora_de_controle, fill = factor(lambda))) +
    geom_boxplot() +
    labs(
      x = "Valores de Φ",
      y = "Fração de pontos fora de controle",
      fill = "Valores de λ",
        title = "Fração de pontos fora de controle por Φ e λ"
    ) +
    theme.base
  ) %>%
  plotly.base %>%
  layout(boxmode = 'group')

Carta CUSUM

cumsum_qcc <- function(amostra_inicial_, se.shift, amostra_subsequente_ = NULL) {
  intervalo_de_decisao <- 5
  
  arguments <- list(
    data = amostra_inicial_,
    se.shift = se.shift,
    decision.interval = intervalo_de_decisao,
    plot = F
  )
  
  if (!is.null(amostra_subsequente_)) {
    arguments$newdata <- amostra_subsequente_
  }
  
  cu <- do.call(cusum, arguments)
  
  registros <- if (is.null(amostra_subsequente_)) {
    seq(1, nH0)
  } else {
    seq(nH0 + 1, nH0 + nH1)
  }
  pos <- cu[["pos"]][registros]
  neg <- cu[["neg"]][registros]
  fora_de_controle_pos = pos > intervalo_de_decisao
  fora_de_controle_neg = neg < -intervalo_de_decisao
  total_fora_de_controle = sum(fora_de_controle_pos) + sum(fora_de_controle_neg)
  fracao_fora_de_controle = total_fora_de_controle / length(registros)
  list(
    pos = pos,
    neg = neg,
    fora_de_controle_pos = fora_de_controle_pos,
    fora_de_controle_neg = fora_de_controle_neg,
    total_fora_de_controle = total_fora_de_controle,
    fracao_fora_de_controle = fracao_fora_de_controle
  )
}

Estimando o valor ótimo de se.shift

Queremos que nas 100 amostras iniciais a probabilidade de um ponto fora de controle seja de 0%. Portanto,

se_shift_otimo_optimize <- function(amostra_inicial) {
  for (x in seq(0, 3, by = 0.1)) {
    if (cumsum_qcc(amostra_inicial, se.shift = x)$total_fora_de_controle == 0) {
      return(x)
    }
  }
}

se_shift_otimo <- cache_dados(
  "cache_se_shift_otimo.RData",
  \() {
    data.frame(
      id = 1:100
    ) %>%
    rowwise() %>%
    mutate(
      amostra_controle = list(barma.sim(
        n = nH0,
        phi = phi_parametro,
        seed = id
      )),
      ultima_observacao = amostra_controle[nH0],
      phi_estimado = list(barma.phi_estimado(amostra_controle)),
      residuos_controle = list(barma.residuos(amostra_controle, phi_estimado)),
      minimo = se_shift_otimo_optimize(residuos_controle)
    )
  }
)

ggplotly(
  se_shift_otimo %>%
    ggplot(
      aes(x = 0, y = minimo, group = 0, fill = "red")
    ) +
    geom_boxplot() +
    labs(
      x = element_blank(),
      y = "Valor ótimo de se.shift",
      title = "Valor ótimo de se.shift por simulação"
    ) +
    theme.base +
    theme.no_legend
) %>%
  plotly.base

valor_otimo_se_shift <- mean(se_shift_otimo$minimo)

kable(
  se_shift_otimo %>%
    group_by() %>%
    summarise(
      "Média" = mean(minimo),
      "Desvio padrão" = sd(minimo)
    ),
  caption = "Valor ótimo de se.shift",
  align = 'c',
  digits = 3
)

Valor ótimo de se.shift

Média	Desvio padrão
0.777	0.29

Simulação de amostras subsequentes

cumsum_amostras_subsequentes_gerador <- function() {
  amostra_controle <- barma.sim(
    n = nH0,
    phi = phi_parametro,
    seed = 1337
  )
  
  ultima_observacao <- amostra_controle[nH0]
  phi_estimado <- barma.phi_estimado(amostra_controle)
  residuos_controle <- barma.residuos(amostra_controle, phi_estimado)
  
  data.frame(
      novo_phi = seq(0, 0.8, by = 0.1)
    ) %>%
    mutate(
      id = row_number(),
    ) %>%
    rowwise() %>%
    mutate(
      observacao = list(1:nH1),
      dados = list(barma.sim(
        n = nH1,
        phi = novo_phi,
        seed = 1337,
        y.start = ultima_observacao
      )),
      residuos = list(
        barma.residuos(dados, phi_estimado)
      ),
      qcc = list(cumsum_qcc(residuos_controle, valor_otimo_se_shift, residuos)),
      pos = list(qcc$pos),
      neg = list(qcc$neg),
      fora_de_controle_pos = list(qcc$fora_de_controle_pos),
      fora_de_controle_neg = list(qcc$fora_de_controle_neg),
      total_fora_de_controle = qcc$total_fora_de_controle,
      fracao_fora_de_controle = qcc$fracao_fora_de_controle
    )
}


cumsum_amostras_subsequentes <- so_quero_que_funcione(cumsum_amostras_subsequentes_gerador, F)

cumsum_amostras_subsequentes %>%
  group_by(novo_phi) %>%
  summarise(
    mean = mean(fracao_fora_de_controle),
    .groups = "drop"
  )

Cartas

cusum_1 <- function(){
  plot <- cumsum_amostras_subsequentes %>%
    filter(novo_phi <= 0.5) %>%
    unnest(
      cols = c(novo_phi, observacao, pos, neg, residuos, fora_de_controle_pos, fora_de_controle_neg)
    ) %>%
    ggplot() +
    geom_hline(aes(yintercept = 5, color = "Limite de decisão"), linetype = "dashed") +
    geom_hline(aes(yintercept = -5, color = "Limite de decisão"), linetype = "dashed") +
    geom_line(aes(x = observacao, y = pos, color = "N+")) +
    geom_line(aes(x = observacao, y = neg, color = "N-")) +
    geom_point(data = . %>% filter(fora_de_controle_pos),
               aes(x = observacao, y = pos, color = "Fora de controle"), size = 1.5, shape = 4) +
    geom_point(data = . %>% filter(fora_de_controle_neg),
               aes(x = observacao, y = neg, color = "Fora de controle"), size = 1.5, shape = 4) +
    labs(x = "Observação", y = "Caracterísitca", color = "Legenda") +
    scale_color_manual(
      values = c(
        "N+" = adjustcolor("blue", alpha.f = 0.6),
        "N-" = adjustcolor("darkgreen", alpha.f = 0.6),
        "Fora de controle" = adjustcolor("red", alpha.f = 0.3),
        "Fora de controle" = adjustcolor("red", alpha.f = 0.4),
        "Limite de decisão" = adjustcolor("red", alpha.f = 0.5)
      )
    ) +
    labs(
      title = "CUMSUM para resíduos βAR(1)",
    ) +
    theme.base +
    theme(
      legend.position = "bottom",
      legend.title = element_blank()
    ) +
    facet_wrap(~novo_phi)
  plot
}

cusum_2 <- function(){
  plot <- cumsum_amostras_subsequentes %>%
    filter(novo_phi > 0.5) %>%
    unnest(
      cols = c(novo_phi, observacao, pos, neg, residuos)
    ) %>%
    ggplot() +
    geom_hline(aes(yintercept = 5, color = "Limite de decisão"), linetype = "dashed") +
    geom_hline(aes(yintercept = -5, color = "Limite de decisão"), linetype = "dashed") +
    geom_line(aes(x = observacao, y = pos, color = "N+")) +
    geom_line(aes(x = observacao, y = neg, color = "N-")) +
    labs(x = "Observação", y = "Caracterísitca", color = "Legenda") +
    scale_color_manual(
      values = c(
        "N+" = adjustcolor("blue", alpha.f = 0.6),
        "N-" = adjustcolor("darkgreen", alpha.f = 0.6),
        "Limite de decisão" = adjustcolor("red", alpha.f = 0.5)
      )
    ) +
    labs(
      title = "CUMSUM para resíduos βAR(1)",
    ) +
    theme.base +
    theme(
      legend.position = "bottom",
      legend.title = element_blank()
    ) +
    facet_wrap(~novo_phi)
  plot
}

ggp_1 <- ggplotly(cusum_1()) %>%
  plotly.base

ggp_2 <- ggplotly(cusum_2()) %>%
  plotly.base

subplot(ggp_1, ggp_2, nrows = 2, margin = 0.05, heights = c(0.7, 0.3))

Simulação de Monte Carlo

cumsum_monte_carlo_gerador <- function(){
  numero_de_execucoes <- 100
  
  expand.grid(
    k = 1:numero_de_execucoes,
    novo_phi = seq(0.2, 0.6, by = 0.1),
    se_shift = seq(0.1, 1.0, by = 0.1)
  ) %>%
  mutate(
    id = row_number()
  ) %>%
  rowwise() %>%
  mutate(
    amostra_controle = list(barma.sim(
      n = nH0,
      phi = phi_parametro,
      seed = id
    )),
    ultima_observacao = amostra_controle[nH0],
    phi_estimado = barma.phi_estimado(amostra_controle),
    residuos_controle = list(
      barma.residuos(amostra_controle, phi_estimado)
    ),
    dados = list(barma.sim(
      n = nH1,
      phi = novo_phi,
      y.start = ultima_observacao,
      seed = id + 1337E4
    )),
    residuos = list(
      barma.residuos(dados, phi_estimado)
    ),
    qcc = list(cumsum_qcc(residuos_controle, se_shift, residuos)),
    pos = list(qcc$pos),
    neg = list(qcc$neg),
    fora_de_controle_pos = list(qcc$fora_de_controle_pos),
    fora_de_controle_neg = list(qcc$fora_de_controle_neg),
    total_fora_de_controle = qcc$total_fora_de_controle,
    fracao_fora_de_controle = qcc$fracao_fora_de_controle
  )
}

cumsum_monte_carlo <- cache_dados(
  "cache_simulacao_cumsum.RData",
  \() so_quero_que_funcione(cumsum_monte_carlo_gerador)
)

Resumo

Em média, diminuir o se.shift aumenta a sensibilidade do CUSUM para detectar mudanças no processo.

cumsum_monte_carlo_resumo <- cumsum_monte_carlo %>%
  group_by(se_shift, novo_phi) %>%
  summarise(
    mean = mean(fracao_fora_de_controle),
    min = min(fracao_fora_de_controle),
    max = max(fracao_fora_de_controle),
    .groups = "drop"
  )

cumsum_monte_carlo_resumo

Gráficos

ggplotly(
  cumsum_monte_carlo_resumo %>%
      ggplot(aes(x = novo_phi, y = mean, color = factor(se_shift))) +
      geom_line() +
      geom_point() +
      labs(
        x = "Valores de Φ",
        y = "Fração de pontos fora de controle",
        color = "Valores de se.shift",
        title = "Fração de pontos fora de controle por Φ e se.shift"
      ) +
      theme.base
) %>%
  plotly.base

Box plot

ggplotly(
  cumsum_monte_carlo %>%
    ggplot(aes(x = factor(novo_phi), y = fracao_fora_de_controle, fill = factor(se_shift))) +
    geom_boxplot() +
    labs(
      x = "Valores de Φ",
      y = "Fração de pontos fora de controle",
      fill = "Valores de se.shift",
        title = "Fração de pontos fora de controle por Φ e se.shift"
    ) +
    theme.base
  ) %>%
  plotly.base %>%
  layout(boxmode = 'group')

EWMA x CUSUM

estatisticas_resumo_combinadas <- bind_rows(
  cumsum_monte_carlo_resumo %>%
    filter(se_shift > 0.7) %>%
    mutate(algoritmo = "CUSUM"),
  ewma_monte_carlo_resumo %>%
    mutate(algoritmo = "EWMA")
)

ggplotly(
  estatisticas_resumo_combinadas %>%
    ggplot() +
    geom_line(
      aes(x = novo_phi, y = mean, color = factor(lambda), linetype = algoritmo),
      data = . %>% filter(algoritmo == "EWMA")
    ) +
    geom_point(
      aes(x = novo_phi, y = mean, color = factor(lambda), shape = algoritmo),
      data = . %>% filter(algoritmo == "EWMA")
    ) +
    geom_line(
      aes(x = novo_phi, y = mean, color = factor(se_shift), linetype = algoritmo),
      data = . %>% filter(algoritmo == "CUSUM")
    ) +
    geom_point(
      aes(x = novo_phi, y = mean, color = factor(se_shift), shape = algoritmo),
      data = . %>% filter(algoritmo == "CUSUM")
    ) +
    labs(
      x = "Valores de Φ",
      y = "Fração de pontos fora de controle",
      linetype = "Algoritmo",
      title = "Fração de pontos fora de controle",
      color = "Valores dos parâmetros",
      shape = "Algoritmo"
    ) +
    theme.base
) %>%
  plotly.base

estatisticas_combinadas <- bind_rows(
  cumsum_monte_carlo %>%
    filter(se_shift > 0.7) %>%
    mutate(algoritmo = "CUSUM"),
  ewma_monte_carlo %>%
    mutate(algoritmo = "EWMA")
)

ggplotly(
  estatisticas_combinadas %>%
    ggplot() +
    geom_boxplot(
      aes(x = factor(novo_phi), y = fracao_fora_de_controle, fill = factor(lambda), shape = algoritmo),
      data = . %>% filter(algoritmo == "EWMA")
    ) +
    geom_boxplot(
      aes(x = factor(novo_phi), y = fracao_fora_de_controle, fill = factor(se_shift), shape = algoritmo),
      data = . %>% filter(algoritmo == "CUSUM")
    ) +
    labs(
      x = "Valores de Φ",
      y = "Fração de pontos fora de controle",
      fill = "Valores dos parâmetros",
      linetype = "Algoritmo",
      title = "Fração de pontos fora de controle",
      color = "Valores dos parâmetros"
    ) +
    facet_wrap(~algoritmo) +
    theme.base
) %>%
  plotly.base %>%
  layout(boxmode = 'group')

E se… combinarmos os dois algoritmos?

comb_monte_carlo_gerador <- function(numero_de_execucoes) {
  numero_de_execucoes <- 100
  
  expand.grid(
    k = 1:numero_de_execucoes,
    novo_phi = seq(0.2, 0.4, by = 0.1),
    se_shift = seq(0.7, 1.1, by = 0.1),
    lambda = seq(0.1, 0.4, by = 0.1)
  ) %>%
  mutate(
    id = row_number(),
    "(se;λ)" = factor(paste0(se_shift, ";", lambda))
  ) %>%
  rowwise() %>%
  mutate(
    amostra_controle = list(barma.sim(n = nH0, phi = phi_parametro, seed = id)),
    ultima_observacao = amostra_controle[nH0],
    phi_estimado = barma.phi_estimado(amostra_controle),
    residuos_controle = list(barma.residuos(amostra_controle, phi_estimado)),
    dados = list(barma.sim(n = nH1, phi = novo_phi, y.start = ultima_observacao, seed = id + 1337E4)),
    residuos = list(barma.residuos(dados, phi_estimado)),
    ewma = list(ewma_qcc(residuos_controle, lambda, residuos)),
    cumsum = list(cumsum_qcc(residuos_controle, se_shift, residuos)),
    fora_de_controle = list(
      ewma$fora_de_controle & (cumsum$fora_de_controle_pos | cumsum$fora_de_controle_neg)
    ),
    fracao_fora_de_controle = sum(fora_de_controle) / nH1
  )
}

comb_monte_carlo <- cache_dados(
  "cache_simulacao_comb.RData",
  \() so_quero_que_funcione(comb_monte_carlo_gerador)
)

Resultados

comb_monte_carlo_resumo <- comb_monte_carlo %>%
  group_by(`(se;λ)`, novo_phi, se_shift, lambda) %>%
  summarise(
    mean = mean(fracao_fora_de_controle),
    min = min(fracao_fora_de_controle),
    max = max(fracao_fora_de_controle),
    .groups = "drop"
  )

comb_monte_carlo_resumo

ggplotly(
  comb_monte_carlo_resumo %>%
    ggplot(aes(x = novo_phi, y = mean, color = `(se;λ)`)) +
    geom_line() +
    geom_point() +
    labs(
      x = "Valores de Φ",
      y = "Fração de pontos fora de controle",
      color = "Valores dos parâmetros (se;λ)",
      title = "Fração de pontos fora de controle por Φ"
    ) +
    theme.base
) %>%
  plotly.base

A melhor combinação encontrada

melhores <- c("1;0.3")

ggplotly(
  comb_monte_carlo_resumo %>%
    filter(`(se;λ)` %in% melhores) %>%
    ggplot(aes(x = novo_phi, y = mean, color = `(se;λ)`)) +
    geom_line() +
    geom_point() +
    labs(
      x = "Valores de Φλ",
      y = "Fração de pontos fora de controle",
      color = "Valores dos parâmetros (se;λ)",
      title = "Fração de pontos fora de controle por Φ"
    ) +
    theme.base
) %>%
  plotly.base

ggplotly(
  comb_monte_carlo %>%
    filter(`(se;λ)` %in% melhores) %>%
    ggplot(aes(x = factor(novo_phi), y = fracao_fora_de_controle, fill = `(se;λ)`)) +
    geom_boxplot() +
    labs(
      x = "Valores de Φ",
      y = "Fração de pontos fora de controle",
      fill = "Valores de `(se;λ)`",
        title = "Fração de pontos fora de controle por Φ e se.shift"
    ) +
    theme.base
  ) %>%
  plotly.base %>%
  layout(boxmode = 'group')